/** ANAYLSIS_12.C
 *  Analyse Macros zum F-Praktikumsversuch 12: Paarspektrometer
 * 
 *  date:   01-Nov-11
 *  author: Ch.Weinsheimer 
 */

#include <iostream>

/**
 * User Defined Functions
 */

/// Define gauss function for fitting
double gaussian(double* vars, double* pars) {
    return (pars[3]+pars[0]*TMath::Gaus(vars[0],pars[1],pars[2]));
}

/// Define linear function for fitting
double lin(double* vars, double* pars) {
    return (pars[0]*vars[0]+pars[1]);
}

/// Define gauss+background function for fitting
double gauss_back(double* vars, double* pars) {
  return pars[0]*TMath::Gaus(vars[0],pars[1],pars[2])+
         pars[3]*TMath::Gaus(vars[0],pars[4],pars[5])+
         pars[6]*vars[0]+pars[7];
}

/// Read data from file and plot a histogram
void histo(char* filename)
{
    int const N = 8192;  

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);

    // Create 1-D Histogram
    TH1F* Hist = gROOT->FindObject("Hist");
    if(!Hist) Hist = new TH1F("Hist", "Histogramm", N, 0, N);
    else Hist->Reset("M");

    // Fill histogramm
    for(Int_t i=0; i < N; i++){
      Hist->Fill(channel[i],value[i]);
    }

    Hist->Rebin(2);

    // Setze Styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);

    // Setze Beschriftungen und Styles
    Hist->SetTitle("MCA Spectrum");
    Hist->GetXaxis()->SetTitle("channel");
    Hist->GetYaxis()->SetTitle("counts");

    // Erzeuge das Fenster
    TCanvas* Canva = gROOT->Find Object("Canva");
    if(!Canva) Canva = new TCanvas("Canva", "Versuch_12", 1024, 768);
    else Canva->Clear();
    Canva->cd(1); // Waehle das Feld an

    // Male Graph
    Hist->Draw("PL");
}



/// Read data from file and plot a histogram
void fitPeaks(char* filename,int const npeaks)
{
    int const N = 8192;  

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);

    // Create 1-D Histogram
    TH1F* fHist = gROOT->FindObject("fHist");
    if(!fHist) fHist = new TH1F("fHist", "Fitted Histogramm", N, 0, N);
    else fHist->Reset("M");

    // Fill histogramm
    for(Int_t i=0; i < N; i++){
      fHist->Fill(channel[i],value[i]);
    }

    fHist->Rebin(2);


    // Erzeuge das Fenster
    gROOT->SetStyle("Plain");
    TCanvas* fCanva = gROOT->Find Object("fCanva");
    if(!fCanva) fCanva = new TCanvas("fCanva", "Fitted Versuch_12", 1024, 768);
    else fCanva->Clear();
    fCanva->cd(1); // Waehle das Feld an

    // Setze Styles
    gStyle->SetOptStat(0);
    //gStyle->SetOptFit(111);

    fHist->Draw("PL");


    // Declare Fit-Parameters and Fit-Funtions
    double mean[npeaks], stdv[npeaks], norm[npeaks], offset[npeaks], min[npeaks], max[npeaks];
    
    TF1 fit[npeaks];
    //    TH1F* hc[npeaks];

    for (unsigned int j = 0; j < npeaks; j++) {

      cout << "Setting Parameters for Peak " << j + 1 << ": " << endl << "******************************" << endl;
      cout << "Position   : "; cin >> mean[j];
      cout << "Lower Bound: "; cin >> min[j];
      cout << "Upper Bound: "; cin >> max[j];
        
      char name[2];
      char histname[4];
      sprintf( name,"%d", j);
      sprintf( histname,"hc%d",j );

      //      hc[j] = (TH1F*)fHist->Clone(histname);

      fit[j] = new TF1(name,gaussian,min[j],max[j],4);
      fit[j].SetParNames("Norm", "Mean", "StdDev", "Offset");
      fit[j].SetParameters(100,mean[j],100,1);
      fit[j].SetLineColor(kRed);

      //      hc[j]->Fit(&fit[j],"R","sames");
      fHist->Fit(&fit[j],"R");
      cout << "Chi-Square: " << fit[j].GetChisquare() << "/" << fit[j].GetNDF() << endl;

    }
    // Canva->Update();



    // Setze Beschriftungen und Styles
    fHist->SetTitle("MCA Spectrum");
    fHist->GetXaxis()->SetTitle("channel");
    fHist->GetYaxis()->SetTitle("counts");
    


    // Male Graph
    fHist->Draw("PL");
    for (j=0; j<npeaks; j++) {
      //hc[j]->DrawClone("Same");
      fit[j].DrawClone("Same");
    }
}


void energyCal(char* filename)
{
    int const N = 6;  

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    Double_t energy[N], channel[N], dchannel[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%lf %lf %lf", &energy[i], &channel[i], &dchannel[i]);
        cout << energy[i] << "   " << channel[i] << "   " << dchannel[i] << endl;
    }
    fclose(file);

    // Erzeuge TGraph repsaesentation
    TGraphErrors* eGraph = gROOT->FindObject("eGraph");
    if(!eGraph) eGraph = new TGraphErrors(N,channel,energy,dchannel,0);
    else Graph->Reset("M");

    // Erzeuge das Fenster
    gROOT->SetStyle("Plain");
    TCanvas* eCanva = gROOT->Find Object("eCanva");
    if(!eCanva) eCanva = new TCanvas("eCanva", "Energy Calibration", 1024, 768);
    else eCanva->Clear();
    eCanva->cd(1); // Waehle das Feld an


    // Setze Styles
    gStyle->SetOptStat(0);
    gStyle->SetOptFit(0);

    TF1 eFunc("eFunc", lin, 0, 1500, 2);
    eFunc.SetParNames("Slope", "Offset");
    eFunc.SetParameters(4,1);
    eFunc.SetLineColor(kRed);
   
    eGraph->Fit(&eFunc,"0");

    // Setze Beschriftungen und Styles
    gStyle->SetOptFit(11);
    eGraph->SetMarkerStyle(2);
    eGraph->SetMarkerSize(2);
    eGraph->SetLineWidth(2);
    eGraph->SetTitle("Energy Calibration");
    eGraph->GetXaxis()->SetTitle("Channel");
    eGraph->GetYaxis()->SetTitle("Energy [keV]");

    eGraph->Draw("APE");
    eFunc.DrawClone("Same");
    
}


/// Plot Time Calibration Measurement
void timecal()
{
  Float_t peaks[6] = {221.0, 269.0, 322.0, 383.0, 422.0, 478.0};
  Float_t dpeaks[6] = {5.0, 5.0, 5.0, 5.0, 5.0, 5.0};
  Float_t times[6] = {10.0, 20.0, 30.0, 40.0, 50.0, 60.0};
  Float_t dtimes[6] = {0.5,0.5,0.5,0.5,0.5,0.5};

  // Create timegraph
  TGraphErrors * gr = new TGraphErrors(6,peaks,times,dpeaks,dtimes);
  gr->SetMarkerColor(kRed);
  gr->SetMarkerStyle(20);
  gr->SetMarkerSize(0.7);

  TF1 tFunc("lin", lin, 200, 500, 2);
  tFunc.SetParNames("Slope", "Offset");
  tFunc.SetParameters(0.2,0);
  tFunc.SetLineColor(kBlue);
  tFunc.SetLineWidth(2);
 
  gr->Fit(&tFunc);



  // Setze Styles
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  // Setze Beschriftungen und Styles
  gr->SetTitle("Time Calibration");
  gr->GetXaxis()->SetTitle("channel");
  gr->GetYaxis()->SetTitle("time / ns");

  // Erzeuge das Fenster
  TCanvas* Canva = gROOT->Find Object("Canva");
  if(!Canva) Canva = new TCanvas("Canva", "Versuch_12 TimeCal", 1024, 768);
  else Canva->Clear();
  Canva->cd(1); // Waehle das Feld an

  // Male Graph
  gr->Draw("APE");
  tFunc.Draw("Same");
}


/// Read data from file and plot a histogram
void polonium(char* filename)
{
    int const N = 8192;
    int const ub = 4200;
    int const lb = 250;

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);

    // Create 1-D Histogram
    TH1F* hp = new TH1F("hp", "Fitted Histogramm", N, 0, N);

    // Fill histogramm
    for(Int_t i=250; i < N; i++){
      hp->Fill(channel[i],value[i]);
    }

    hp->Rebin(2);


    // Erzeuge das Fenster
    gROOT->SetStyle("Plain");
    TCanvas* c1 = new TCanvas("c1", "Fitted Versuch_12", 1024, 768);
    c1->cd(1); // Waehle das Feld an

    // Setze Styles
    gStyle->SetOptStat(0);
    gStyle->SetOptFit(111);


    // Declare Fit-Parameters and Fit-Funtions
    TF1 background("background","[0]*x+[1]",lb,ub);
    background.SetLineWidth(2);
    background.SetLineStyle(2);
    background.SetLineColor(kBlue);

    TF1 signal("signal","[0]*TMath::Gaus(x,[1],[2])+[3]*TMath::Gaus(x,[4],[5])",lb,ub);
    signal.SetLineWidth(2);
    signal.SetLineStyle(2);
    signal.SetLineColor(kRed);

    TF1 fit("fit",gauss_back,lb,ub,8);
    fit.SetLineWidth(2);
    fit.SetLineColor(kRed);
    double norm1 = 800.0; double mean1 = 2250.0; double stddev1 = 100.0;
    double norm2 = 300.0; double mean2 = 3500.0; double stddev2 = 150.0;
    double m = 0.027; double offset = 1.0;
    fit.SetParameters(norm1,mean1,stddev1,norm2,mean2,stddev2,m,offset);
    fit.SetParNames("Norm1","Mean1","StdDev1","Norm2","Mean2","StdDev2","m","Offset");

    // Perform Fit
    hp->Fit(&fit,"R");
    cout << "Chi-Square: " << fit.GetChisquare() << "/" << fit.GetNDF() << endl;


    background.SetParameters( fit.GetParameter(6), fit.GetParameter(7) );
    signal.SetParameters( fit.GetParameter(0), fit.GetParameter(1), fit.GetParameter(2),
			  fit.GetParameter(3), fit.GetParameter(4), fit.GetParameter(5)
			);

    // Setze Beschriftungen und Styles
    hp->SetTitle("MCA Spectrum");
    hp->GetXaxis()->SetTitle("channel");
    hp->GetYaxis()->SetTitle("counts");
    
    hp->GetYaxis()->SetRangeUser(0,1050);

    // Male Graph
    hp->Draw("PL");
    background.DrawClone("Same");
    signal.DrawClone("Same");


}
